SCP - Centrality

1 Carga de librerías

if (!(require(DT)))
  install.packages("DT")
library(DT)

if (!require(readxl))
  install.packages("readxl")
library(readxl)

# Rfast para matriz de varianza combinada
if (!requireNamespace("Rfast"))
  install.packages("Rfast")
library(Rfast)

if (!requireNamespace("ggplot2"))
  install.packages("ggplot2")
library(ggplot2)

if (!requireNamespace("plotly"))
  install.packages("plotly")
library(plotly)

if (!requireNamespace("ggrepel"))
  install.packages("ggrepel")
library(ggrepel)

# GridFCM
library(devtools)
if (!requireNamespace("GridFCM.practicum"))
  install_github("asanfe/GridFCM.practicum", quietly = TRUE)
library(GridFCM.practicum)

# Viridislilte
if (!requireNamespace("viridisLite"))
  install.packages("viridisLite")
library(viridisLite)

# Test para normalidad multivariante
if (!requireNamespace("MVN"))
  install.packages("MVN")
library(MVN)

if (!requireNamespace("ggpattern", quietly = TRUE))
  install.packages("ggpattern")
library(ggpattern)

if (!requireNamespace("factoextra", quietly = TRUE))
  install.packages("factoextra")
library(factoextra)

if (!requireNamespace("cluster", quietly = TRUE))
  install.packages("cluster")
library(cluster)

if (!requireNamespace("RColorBrewer", quietly = TRUE))
  install.packages("RColorBrewer")
library(RColorBrewer)

if (!requireNamespace("rcartocolor", quietly = TRUE))
  install.packages("rcartocolor")
library(rcartocolor)

2 Importación y resumen

2.1 Importación del objeto RDA

# Objetos de sesión de ejemplo de la PEC
path <- '../CentralityTest/data.RData'
load(path)
samples.raw.df <- data

samples.df <- data.frame(ID = samples.raw.df$dataset$ID,
                         gender = samples.raw.df$dataset$gender,
                         age = samples.raw.df$dataset$age,
                         edu = samples.raw.df$dataset$edu,
                         status = "Error")

for(i in 1:nrow(samples.df)) {
  id <- samples.df$ID[i]  
  
  tryCatch({
    wimp <- data$grids[[id]]$WT  # Accede al wimp asociado al ID
    wphm <- GridFCM.practicum::ph_index(wimp = wimp, method = "wnorm", std = "none")

    samples.df$status[i] <- "Ok"
  }, error = function(e){
    cat("Error procesando ID:", id, "\n")
  })
}
## Error procesando ID: P0166567 
## Error procesando ID: P0512115 
## Error procesando ID: P0606903 
## Error procesando ID: P0666512 
## Error procesando ID: P0870223 
## Error procesando ID: P0910424 
## Error procesando ID: P1102149 
## Error procesando ID: P1123165 
## Error procesando ID: P1140623 
## Error procesando ID: P1312902 
## Error procesando ID: P1426704 
## Error procesando ID: P1554581 
## Error procesando ID: P1891931 
## Error procesando ID: P1933446
# Calcula las frecuencias de status de procesamiento
freqs.status <- table(samples.df$status, useNA = "no")

# Calcula los porcentajes
percent.status <- prop.table(freqs.status) * 100

results.summary <- data.frame(
  ResultadoCarga = names(freqs.status),
  Casos = as.integer(freqs.status),
  Porcentaje = round(100 * prop.table(freqs.status), 3)
)

results.summary <- results.summary[, c("ResultadoCarga", "Casos", "Porcentaje.Freq")]

2.2 Resultado de procesamiento

# Mostrar la tabla
DT::datatable(results.summary, options = list(pageLength = 5))

2.3 Resumen de sujetos

DT::datatable(samples.df)

2.4 Wimp del sujeto a presentar

#path <- 'Wimp_Ejemplo.xlsx'
#opr <- TRUE
#wimp <- GridFCM.practicum::importwimp(path = path, opr = opr, sheet = 1)

# Objetos de sesión de ejemplo de la PEC
path <- '../CentralityTest/data.RData'
load(path)
samples.raw.df <- data

# Selección de un sujeto
id.sample <- params$id.sujeto.entrada
wimp <- data$grids[[id.sample]]$WT

bertin(wimp$openrepgrid, colors = c("palegreen", "darkgreen"))

3 Exploración de la Wimp

3.1 Digrafo del Self

# Digraph
GridFCM.practicum::digraph(wimp, layout = "rtcircle")

3.2 Digrafo del Ideal

# Digraph
GridFCM.practicum::idealdigraph(wimp, layout = "rtcircle")

3.3 E/S de los constructos. Método Simple

c.io.test <- GridFCM.practicum::degree_index(wimp, method = "simple")
c.io.test <- c.io.test[, c(2,1,3)]
c.io.test <- cbind(c.io.test, Diff = (c.io.test[, 2] - c.io.test[, 1]))
c.io.test.r <- round(c.io.test, 3)
DT::datatable(c.io.test.r)

3.4 E/S de los constructos. Método wnorm

c.io.test <- GridFCM.practicum::degree_index(wimp, method = "wnorm")
c.io.test <- c.io.test[, c(2,1,3)]
c.io.test <- cbind(c.io.test, Diff = (c.io.test[, 2] - c.io.test[, 1]))
c.io.test.r <- round(c.io.test, 3)
DT::datatable(c.io.test.r)

3.5 Ejemplo de salida de P-H index

test.wphm <- GridFCM.practicum::ph_index(wimp = wimp, method = "wnorm", std = FALSE)
DT::datatable(test.wphm)

4 Distancia de Mahalanobis y distribución de datos

4.1 Test de Mardia para análisis multivariante

Llevamos a cabo previamente un test de Mardia para constrastar la normalidad multivariante de los datos, a fin de determinar la pertinencia del punto de corte basado en adecuación a distribución Chi-cuadrado de distancia de Mahalanobis

# Test de Mardia

test.result <- mvn(data = test.wphm, mvnTest = "mardia")

print(test.result)
## $multivariateNormality
##              Test          Statistic            p value Result
## 1 Mardia Skewness    8.0225087103244 0.0907571427349819    YES
## 2 Mardia Kurtosis -0.110858005728675  0.911728946847989    YES
## 3             MVN               <NA>               <NA>    YES
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling     p        0.4564    0.2183    YES   
## 2 Anderson-Darling     h        0.6781    0.0567    YES   
## 
## $Descriptives
##    n         Mean    Std.Dev     Median        Min       Max        25th
## p 12 3.839090e-01 0.22974792 0.33105454  0.1285649 0.8678129  0.18909750
## h 12 4.626635e-18 0.09769977 0.02249885 -0.2003469 0.1628488 -0.01151727
##         75th       Skew   Kurtosis
## p 0.56488189  0.5943136 -0.9142834
## h 0.04258711 -0.6357197 -0.2658073

4.2 Test de resultado de la función

test.wmahalanobis <- GridFCM.practicum::mahalanobis_index(wimp = wimp, method = "wnorm", std = FALSE, sign.level = 0.2)
DT::datatable(test.wmahalanobis)

4.3 Distribución Chi-cuadrado

# Definimos los grados de libertad para la distribución chi-cuadrado
df <- 2

# Generamos los valores de la distribución
x <- seq(qchisq(0.001, df), qchisq(0.999, df), length.out = 1000)
y <- dchisq(x, df)

# Calculamos los puntos de corte para el 20% superior
sign.level <- 0.2
cut_high <- qchisq(1- sign.level, df)

# Dataframe para la gráfica
data <- data.frame(x = x, y = y)

ggplot(data, aes(x = x, y = y)) + 
  geom_line() + 
  geom_ribbon(data = data %>% filter(x > cut_high), aes(ymax = y), ymin = 0, fill = 'salmon', alpha = 0.5) +
  geom_vline(xintercept = cut_high, color = "red", linetype = "dashed") +
  labs(title = 'Distribución Chi-cuadrado con puntos de corte del 80%', x = 'Valor', y = 'Densidad') +
  theme_minimal()

4.4 Gráfica de barras de distancias de Mahalanobis y punto de corte

# Distancia de Mahalanobis
test.bp.wmahalanobis <- GridFCM.practicum::mahalanobis_index(wimp = wimp, method = "wnorm", std = FALSE)
test.wmahalanobis.df <- as.data.frame(test.bp.wmahalanobis)

# Colores de los constructos


#test.wmahalanobis.df$constructo <- rownames(test.wmahalanobis)
test.wmahalanobis.df$constructo <- wimp$constructs$right.poles

# Valoración del ideal
test.wmahalanobis.df$idealdirect <- wimp$ideal$direct

# Columna para identificar constructos dilemáticos
#test.wmahalanobis.df$fill.color <- ifelse(test.wmahalanobis.df$idealdirect == 4, "yellow2", "honeydew")
test.wmahalanobis.df$fill.color <- construct_colors(wimp= wimp, mode = "red/green")

# Ordenamos las barras en orden decreciente
test.wmahalanobis.df <- test.wmahalanobis.df %>%
  arrange(desc(m.dist))

# Convertimos 'constructo' en un factor con los niveles en el orden deseado
test.wmahalanobis.df$constructo <- factor(test.wmahalanobis.df$constructo, levels = test.wmahalanobis.df$constructo)

# Punto de corte distribución Chi-Cuadrado
sign.level <- 0.2
df <- ncol(test.wphm)
chi.square.cutoff <- qchisq(1 - sign.level, df)
#media_m_dist <- mean(test.wmahalanobis.df$m.dist)

4.5 Constructos supraordenados

# Crear el histograma de constructos supraordenados

# Filtramos por los constructos donde el valor de 'h' es mayor que cero
test.wmahalanobis.df.sup <- test.wmahalanobis.df %>%
  filter(h > 0)

bar_plot <- ggplot(test.wmahalanobis.df.sup, aes(x = constructo, y = m.dist, fill = fill.color)) +
  geom_bar(stat = "identity", color = "black", linewidth = 0.25) +
  geom_hline(yintercept = chi.square.cutoff, linetype = "dashed", color = "darkgreen", linewidth = 1) +
  scale_fill_identity() + # Usa los colores asignados directamente
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none" # Ocultar leyenda para fill
  ) +
  labs(x = "Constructos con h > 0", y = "Distancia de Mahalanobis", title = "Constructos con h>0 por distancia de Mahalanobis")

# Mostramos el gráfico
print(bar_plot)

4.6 Constructos subordinados

# Crear el histograma de constructos subordinados

# Filtramos por los constructos donde el valor de 'h' es menor que cero
test.wmahalanobis.df.sub <- test.wmahalanobis.df %>%
  filter(h < 0)

bar_plot <- ggplot(test.wmahalanobis.df.sub, aes(x = constructo, y = m.dist, fill = fill.color)) +
  geom_bar(stat = "identity", color = "black", linewidth = 0.25) +
  geom_hline(yintercept = chi.square.cutoff, linetype = "dashed", color = "darkgreen", linewidth = 1) +
  scale_fill_identity() + # Usa los colores asignados directamente
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none" # Ocultar leyenda para fill
  ) +
  labs(x = "Constructos con h < 0", y = "Distancia de Mahalanobis", title = "Constructos con h<0 por distancia de Mahalanobis")

# Mostramos el gráfico
print(bar_plot)

4.7 Distribución de valores en P

4.7.1 Distribución de valores de P

# Crear el histograma de valores en P

test.wmahalanobis.df.sortP <- test.wmahalanobis.df %>% 
  arrange(desc(p))

mean.p <- mean(abs(test.wmahalanobis.df.sortP$p))

# Convertir 'constructo' en un factor para mantener el orden en el gráfico
test.wmahalanobis.df.sortP$constructo <- factor(test.wmahalanobis.df.sortP$constructo, 
                                          levels = test.wmahalanobis.df.sortP$constructo)

bar_plot <- ggplot(test.wmahalanobis.df.sortP, aes(x = constructo, y = p, fill = fill.color)) +
  geom_bar(stat = "identity", color = "black", linewidth = 0.25) +
  geom_hline(yintercept = mean.p, linetype = "dashed", color = "darkgreen", linewidth = 1) +
  scale_fill_identity() + # Usa los colores asignados directamente
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none" # Ocultar leyenda para fill
  ) +
  labs(x = "Constructos", y = "Valor de P", title = "Constructos por distancia en P")

# Mostramos el gráfico
print(bar_plot)

4.7.2 Test Shapiro-Wilk (muestras pequeñas o moderadas) para variable P

# Test de normalidad de Saphiro-Wilk
norm.test <- shapiro.test(test.wmahalanobis.df$p)

# Imprime el resultado
print(norm.test)
## 
##  Shapiro-Wilk normality test
## 
## data:  test.wmahalanobis.df$p
## W = 0.90483, p-value = 0.1831
p.value <- 0.05
norm.test.result <- norm.test$p.value > p.value 

De acuerdo con el resultado de la prueba, la normalidad de la distribución de los valores de P es: TRUE

4.7.3 Gráfica Cuantil-Cuantil

datos <- test.wmahalanobis.df$p

# Gráfica q-q para comprobar la normalidad
qqnorm(datos)
qqline(datos, col = "red")

4.7.4 Punto de corte basado en distribución normal de P

# Media y desviación típica de la distribución
mean.p <- mean(test.wmahalanobis.df$p)
sd.p <- sd(test.wmahalanobis.df$p)

# Definir el rango de valores para X basado en la media y desviación típica
x.values <- seq(from = mean.p - 4 * sd.p, to = mean.p + 4 * sd.p, length.out = 1000)

# Crear un dataframe con los valores de X y la densidad de una distribución normal para esos valores
norm.df <- data.frame(x = x.values, y = dnorm(x.values, mean = mean.p, sd = sd.p))

# Punto de corte
cut.low <- qnorm(0.15, mean = mean.p, sd = sd.p)

plot <- ggplot(norm.df, aes(x = x, y = y)) +
  geom_line() +
  geom_vline(xintercept = cut.low, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = mean.p, linetype = "dashed", color = "blue", linewidth = 1) +
  geom_area(data = subset(norm.df, x <= cut.low), fill = "lightblue", alpha = 0.2) +
  theme_bw() +
  theme(
    panel.grid.major = element_line(linewidth = 0.5, linetype = 'solid', colour = "lightgrey"),
    panel.grid.minor = element_blank(),
    legend.position = "none"
  ) +
  scale_x_continuous(name = "Valor de P") +
  scale_y_continuous(name = "Densidad") +
  labs(title = paste("Distribución Normal con Media en", round(mean.p, 2),
                     "y Punto de Corte en", round(cut.low, 2)))
# Mostrar la gráfica
print(plot)

5 Representación en espacio P-H

5.1 Sin estandarización - plotly sin marcar área no viable ni constructos centrales

wp1.grph <- GridFCM.practicum::graph_ph(wimp = wimp, method = "wnorm", std = 'none', sign.level = 0.2,
                                        mark.nva = FALSE, mark.cnt = FALSE)
wp1.grph

5.2 Espacio PH con coloreado de área no útil y marcado de outliers. Función de representación

5.2.1 Sin estandarización

wp1.grph <- GridFCM.practicum::graph_ph(wimp = wimp, method = "wnorm", std = 'none', sign.level = 0.2,
                                        mark.nva = TRUE, mark.cnt = TRUE, show.points = TRUE)
wp1.grph

5.2.2 Con estandarización basada en aristas

wp1.grph <- GridFCM.practicum::graph_ph(wimp = wimp, method = "wnorm", std = 'edges', sign.level = 0.2,
                                        mark.nva = TRUE, mark.cnt = TRUE, show.points = TRUE)
wp1.grph

5.2.3 Sin estandarización, marcando área de coordenadas no viables. Sin puntos

wp1.grph <- GridFCM.practicum::graph_ph(wimp = wimp, method = "wnorm", std = "none", sign.level = 0.2,
                                        mark.nva = TRUE, mark.cnt = TRUE, show.points = FALSE)
wp1.grph
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

5.2.4 Sin estandarización, sin marcar área de coordenadas no viables. Sin puntos

wp1.grph <- GridFCM.practicum::graph_ph(wimp = wimp, method = "wnorm", std = "none", sign.level = 0.2,
                                        mark.nva = FALSE, mark.cnt = TRUE)

wp1.grph

6 Otros métodos de centralidad

6.1 Primer eigenvectorsobre matriz de implicaciones

eigen_index <- function(wimp){

  # Matriz de adyacencia
  adj.matrix <- wimp$scores$implications

  # Vectores y valores propios
  results <- eigen(adj.matrix)

  # Extrae el primer vector propio, asociado a la varianza máxima de los datos
  eigenvector.principal <- results$vectors[,1]

  # Creamos un marco de datos con los nombres de los constructos y las componentes del primer eigenvector
  df.centrality <- data.frame(
    constructs = wimp$constructs$constructs,
    leftpoles = wimp$constructs$left.poles,
    rightpoles = wimp$constructs$right.poles,
    firsteigenvector = eigenvector.principal
  )

  return(df.centrality)
}
cent.evalues.df <- eigen_index(wimp)
DT::datatable(cent.evalues.df)

6.2 Análisis de componentes principales (PCA)

# Foco del PCA
adj.matrix <- wimp$scores$implications

# Análisis de componentes principales
pca.result <- prcomp(adj.matrix, center = TRUE, scale = TRUE)

# Extraer los dos primeros componentes principales
pca.comp <- as.data.frame(pca.result$x[, 1:2])
pca.comp$constructs <- wimp$constructs$constructs

# Crea la gráfica de dispersión
pca.plot <- plot_ly(data = pca.comp, x = ~PC1, y = ~PC2, type = 'scatter', mode = 'markers',
                    hoverinfo = 'text+x+y',
                    marker = list(size = 10, opacity = 0.8)) %>%
            layout(title = 'PCA de matriz de adyacencia',
                   xaxis = list(title = 'PCA 1'),
                   yaxis = list(title = 'PCA 2'),
                   hovermode = 'closest',
                   plot_bgcolor = "white",
                   font = list(family = "Arial"),
                   showlegend = FALSE) %>%
            # Add annotations for each point
            add_annotations(data = pca.comp, x = ~PC1, y = ~PC2, text = ~constructs,
                            showarrow = FALSE, xanchor = 'center', yanchor = 'bottom', font = list(size = 12))

# Muestra la gráfica
pca.plot

7 Análisis por conglomerados

7.1 Determinación de número óptimo de conglomerados

k <- test_optimal_num_clusters(wimp = wimp, method = "wnorm", std = "none")
k
## [1] 3

Tenemos un número máximo de 3 conglomerados en nuestros datos.

7.2 Representación de números de conglomerados óptimo

Adecuación de cohesión y separación de cada punto según pertenezca a distintos conglomerados:

# Lista que albergará las distintas gráficas de silueta
lista.graf.sil <- list()

# Valores intermedios que ya calcula .optimal.num.clusters
max.clusters <- length(wimp$constructs$constructs) - 1
ph.mat <- GridFCM.practicum::ph_index(wimp = wimp, method = "wnorm", std = "none")
#rownames(test.dist) <- as.character(wimp$constructs$right.poles)
#distancias <- dist(test.dist, method = "euclidean")

#------------------------------
# Matriz de disimilaridad modelada como matriz de distancias de Mahalanobis
# Matriz de covarianzas
cov.matrix <- cov(ph.mat)  # Calcula la matriz de covarianza
# Vector de medias
means.vector <- colMeans(ph.mat)

# Inicializa una matriz para guardar las distancias de Mahalanobis
n <- nrow(ph.mat)
dist.mat <- matrix(NA, n, n)

# Calcula la distancia de Mahalanobis entre cada par de filas en ph.mat
for (i in 1:n) {
  for (j in i:n) {
    diff <- ph.mat[i, ] - ph.mat[j, ]
    dist.mat[i, j] <- sqrt(t(diff) %*% solve(cov.matrix) %*% diff)
    dist.mat[j, i] <- dist.mat[i, j]  # La matriz es simétrica
  }
}

# Hacemos 0 en la diagonal
diag(dist.mat) <- 0

row.names(dist.mat) <- row.names(ph.mat)
colnames(dist.mat) <- row.names(ph.mat)

#---------------------------

# Preparamos una lista con diversas representaciones gráficas de siluetas (de 2 a 10 clústeres)
for(j in 2:min(13, max.clusters)){
  it.pam <- cluster::pam(dist.mat, j, diss = TRUE)
  p <- factoextra::fviz_silhouette(it.pam, label = FALSE, print.summary = FALSE)
  lista.graf.sil[[j-1]] <- p
}

# Organizar los gráficos en una matriz de 4x3, y los presentamos
gridExtra::grid.arrange(grobs = lista.graf.sil, ncol = 3, nrow = 4)

7.2.1 Dendrograma

#Dendrograma
dendron.plot <- constructs_dendrogram(wimp = wimp)
## Registered S3 method overwritten by 'dendextend':
##   method       from   
##   text.pvclust pvclust
print(dendron.plot)

7.2.2 ClusPlot

act.cex <- par("cex")
par(cex = 0.8)

# Calculamos el objeto de partición PAM para los k clústeres obtenidos anteriormente
opt.pam <- cluster::pam(dist.mat, k, diss = TRUE)

# Colores de los clusters
clus.colors <- carto_pal(n = k, "Peach")

cluster::clusplot(x = dist.mat,
                  clus = opt.pam$clustering,
                  shade = TRUE,
                  color = TRUE,
                  col.clus = clus.colors,
                  col.p = 'darkblue',
                  diss = TRUE,
                  labels = 3)